#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wed Jun  7 15:29:11 2023

@author: jcfq2
"""

import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import numpy as np
import shapely.geometry as sgeom
from cartopy.feature import ShapelyFeature
from scipy.io import readsav
import cartopy.crs as ccrs


import sys
import os
sys.path.insert(0, '../../planetary_values')

from grille_vogt_localtime_dawndusk_all import load_vogt_slices
from grille_vogt_localtime_dawndusk_all import vogt_ofl
from grille_vogt_localtime_dawndusk_all import fit_fieldline_fixed
from grille_vogt_localtime_dawndusk_all import fit_fieldline_scattered

from shift_magnetic import to_magnetic
from shift_magnetic import from_magnetic
from planet_magfield import Bmag




sys.path.insert(0, '../../tss')
from basic.polygonmask import maskeqvalue
from basic.polygonmask import masksubarray
from basic.interpolate_irregular import data_atlatlong
 
sys.path.insert(0, '../../planetary_values')

# In[0]
# set up values

bbox_props1= dict(boxstyle="round,pad=0.15", fc="aliceblue", ec="gray", lw=1)
lt_text=['08h','10h','12h','14h','16h']




def plot_labels_on(ax,white=True,ilabel="",inum=99):
        from matplotlib.patches import Circle
        x0=0
        y0=0
        radii=15
        patches = []
        circle = Circle((x0, y0), radii,color='w',ec='k',zorder=5)
 
        
        if white: 
            ax.add_patch(circle)

            col='w'
            col2='k'
        else:
            col='k'
            col2='w'
        
        txt1=ax.text(s="'Dark' region",x=50, y=80, fontsize=10, fontweight='bold' ,fontstyle='italic',c='w',horizontalalignment='center', verticalalignment='center',zorder=4)
        txt2=ax.text(s="Dusk arcs",x=0, y=-70, fontsize=10, fontweight='bold' ,fontstyle='italic',c='w',horizontalalignment='center', verticalalignment='center',zorder=4)
        txt3=ax.text(s="$\Delta$I",x=-65, y=5, fontsize=10, fontweight='bold' ,fontstyle='italic',c='w',horizontalalignment='center', verticalalignment='center',zorder=4)
        # txt4=ax.text(s="'Active'?",x=-60, y=10, fontsize=10, fontweight='bold' ,fontstyle='italic',c='w',horizontalalignment='center', verticalalignment='center')
        
        if ilabel == "":
            pass
        else:
            txt5=ax.text(s=ilabel,x=-60, y=-135, fontsize=10, fontweight='bold' ,fontstyle='italic',c='k',horizontalalignment='center', verticalalignment='center')

        xdr=np.array([0,-25,-50,-40,-20,25,70])
        ydr=np.array([0,6,50,75,100,125,130])
            
        from scipy import interpolate
        
        xsize=20
        
        tck, u = interpolate.splprep([xdr, ydr], s=0)
        xdr_f,ydr_f = interpolate.splev(np.linspace(0, 1, xsize), tck)
        
        # former boundary of the P region

        # xdr=np.array([0,-25,-50,-60,-75,-100])
        # ydr=np.array([0,6,17,20,23,28])
         
        
        # tck, u = interpolate.splprep([xdr, ydr], s=0)
        # xdr_f2,ydr_f2 = interpolate.splev(np.linspace(0, 1, xsize), tck)
        
        xdr=np.array([0,-25,-50,-75,-100])
        ydr=np.array([0,-18,-26,-30,-35])
         
        
        tck, u = interpolate.splprep([xdr, ydr], s=0)
        xdr_f3,ydr_f3 = interpolate.splev(np.linspace(0, 1, xsize), tck)
        
        
        # print(xdr)
        # print(ydr_f)
        
        ax.plot(xdr_f, ydr_f,c=col,linestyle='-.')
        # ax.plot(xdr_f2, ydr_f2,c=col,linestyle='-.')
        ax.plot(xdr_f3, ydr_f3,c=col,linestyle='-.')
        
        
        txt1.set_bbox(dict(facecolor=col2, alpha=0.2, edgecolor=col2, boxstyle='round,pad=0.2'))
        txt2.set_bbox(dict(facecolor=col2, alpha=0.3, edgecolor=col2, boxstyle='round,pad=0.2'))
        txt3.set_bbox(dict(facecolor=col2, alpha=0.3, edgecolor=col2, boxstyle='round,pad=0.2'))



        if inum < 99:
            xsize=50

            if inum==5:
                xdr2=np.array([150,50,0,-25,-25,-15,0,25,50,75,100,125,150])
                ydr2=np.array([-150,-125,-100,-75,-40,-25,0,28,35,38,35,25,10])

            if inum==4:
                xdr2=np.array([180,80,40,0,0,0,0,25,50,75,100,125,150])
                ydr2=np.array([-150,-125,-100,-75,-40,-25,0,50,60,70,70,52,48])

            if inum==3:
                xdr2=np.array([200,100,50,35,25,15,0,-8,-10,2,45,75,150])
                ydr2=np.array([-150,-110,-103,-75,-40,-25,0,25,50,75,95,95,80])

            if inum==2:
               xdr2=np.array([150,125,115,75,50,30,0,-20,-25,-24,0,50,100,150])
               ydr2=np.array([-120,-110,-100,-70,-40,-25,0,25,50,75,100,125,130,130])

            if inum==1:
               xdr2=np.array([150,125,100,75,50,25,0,-40,-50,-40,-24,0,50,100,150])
               ydr2=np.array([-60,-50,-40,-22,-15,-3,0,25,50,75,100,120,130,135,140])
           
            print('INUM:',inum,xdr2)

            tck, u = interpolate.splprep([xdr2, ydr2], s=1)
            xdr_fA,ydr_fA = interpolate.splev(np.linspace(0, 1, xsize), tck)

            # ax.plot(xdr_fA, ydr_fA,c='k',linestyle='-.')
            ax.fill(xdr_fA, ydr_fA,c='w',zorder=3)

maxint=1.0
minint=0.2
mincount=0
maxcount=3600*100.#3.*3600#hours - many hours due to bleeding of individual pixels across multiple bins
mindate=10
maxdate=30


minrange=0.7
maxrange=1.9
minrange=0.0003
maxrange=0.0018;0.0001
xxrange=6
startx=1

polelat=75
polelong=180


# In[1]
# load data


results='/Users/jcfq2/OneDrive - Northumbria University - Production Azure AD/results/'
planetary='/Users/jcfq2/OneDrive - Northumbria University - Production Azure AD/python/planetary_values/'


idl_file = readsav(results+'jupimage_phase_total_crazytown.idl')
# idl_file = readsav(results+'jupimage_phase_total_crazytown_avgclicks.idl')
# idl_file = readsav(results+'jupimage_phase_total_crazytown_avgclicks.idl')
a=idl_file['a']
a2=idl_file['a2']
b=idl_file['b']
b2=idl_file['b2']

aa=a
aa2=a2
bb=b
bb2=b2


ab=a/b
ab=np.nan_to_num(ab)

   # In[3] set up vogt


# for cml in range(50,329,10):

               


# In[4] 

totalmap=np.sum(ab,axis=0)
# plt.imshow(totalmap)
# plt.show()

totalmapN=totalmap
totalmapS=totalmap
totalmapN[0:70,:]=0
totalmapS[90:,:]=0

totalmap=totalmap/np.mean(totalmap)
totalmapN=totalmapN/np.mean(totalmapN)
totalmapS=totalmapS/np.mean(totalmapS)

lon = np.linspace(360, 0, 360)
lat = np.linspace(-90, 90, 181)
lon2d, lat2d = np.meshgrid(lon, lat)





lt_scaling=[0.8828533281700559, 1.2033086368036778, 1.310723368138646, 1.2050975221749605, 0.8864310989126216]

# using the body of the planet emission for lt
lt_scaling=[0.9573501789864948, 1.2452566354398804, 1.3159839931200386, 1.1695322520269686, 0.8059014121606716]
sys.path.insert(0, '../../planetary_values')



sector1_lon,sector1_lat,sector1_loc,sector1_dis, sector2_lon,sector2_lat,sector2_loc,sector2_dis, sector3_lon,sector3_lat,sector3_loc,sector3_dis, sector4_lon,sector4_lat,sector4_loc,sector4_dis, sector5_lon,sector5_lat,sector5_loc,sector5_dis = load_vogt_slices()


   
# In[4]
sys.path.insert(0, '../../planetary_values')
   
ofl_lat,ofl_lon=vogt_ofl()
    
subsolarlat=0
subsolarlongitude=180#240.3#+180


s_ssl=str(int(subsolarlongitude/10+1)*10)

# print(s_ssl)

sys.path.insert(0, '../../planetary_values')

v = open(os.path.join(sys.path[0],'fieldline_tracing/fieldline_tracing_mapping_contours_kk2009ext_jrm09int_north_sslong'+s_ssl+'.txt'))
north_plotting = np.zeros([28,36,4])

header = v.readline()

for z in range(28):
    for y in range(36):
        line = v.readline()
        north_plotting[z,y,0]=line[0:12]
        north_plotting[z,y,1]=line[14:20]
        null1=line[21:29]
        if any(chr.isdigit() for chr in null1):
            north_plotting[z,y,2]=null1
        else:
            north_plotting[z,y,2]=np.nan
        null2=line[30:38]
        if any(chr.isdigit() for chr in null2):
            north_plotting[z,y,3]=null2
        else:
            north_plotting[z,y,3]=np.nan
                                  
        
    
    # In[5]

    
brightness=np.zeros([xxrange,36])
brightness_diff=np.zeros([xxrange,36])
brightness_diff2=np.zeros([xxrange,36])
    

    
for daf in range(2):

    if daf == 0:
        doallfive=True
    else:
        doallfive=False
    
    
    # for i in range(startx,xxrange+2):
    for i in range(startx,xxrange):
    # for i in range(startx+2,startx+3):
    
        
    
        # cml=360-(75+i*30+15)    
        cml=360-(90+i*20+10)   
        
        
        subsolarlat=0
        subsolarlongitude=cml#240.3#+180
        
        
        s_ssl=str(int(subsolarlongitude/10+1)*10)
        
        # print(s_ssl)
        
        
        v = open(planetary+'fieldline_tracing/fieldline_tracing_mapping_contours_kk2009ext_jrm09int_north_sslong'+s_ssl+'.txt')
        north_plotting = np.zeros([28,36,4])
        
        header = v.readline()
        
        for z in range(28):
            for y in range(36):
                line = v.readline()
                north_plotting[z,y,0]=line[0:12]
                north_plotting[z,y,1]=line[14:20]
                null1=line[21:29]
                if any(chr.isdigit() for chr in null1):
                    north_plotting[z,y,2]=null1
                else:
                    north_plotting[z,y,2]=np.nan
                null2=line[30:38]
                if any(chr.isdigit() for chr in null2):
                    north_plotting[z,y,3]=null2
                else:
                    north_plotting[z,y,3]=np.nan
    
    
        # ax1 = plt.subplot(1, 7, 8-i,projection=ccrs.Orthographic(360-polelong, polelat))
        # ax1 = plt.subplot(3, 6, 6-i,projection=ccrs.Orthographic(360-polelong, polelat))
        mapfull=ab
    
        # mapmap=np.sum(mapfull[90+i*20:90+(i+1)*20,:,:],axis=0)
        mapmap=np.sum(mapfull[75+i*30:75+(i+1)*30,:,:],axis=0)
        mapmapT=np.sum(mapfull[75+1*30:75+(5+1)*30,:,:],axis=0)
        # mapmap=mapmap/np.mean(np.nonzero(mapmapT))
    
        countmap=np.sum(b[75+i*30:75+(i+1)*30,:,:],axis=0)
        totalmap3=np.sum(a[75+i*30:75+(i+1)*30,:,:],axis=0)
        
        
        ax0=1120.572303190426
        ax1=-43.641738021841704#g.parameters[0]
    
        B1=Bmag(lon2d, lat2d,radial=False)
        
        scaleint=B1*ax1+ax0
        scaleint=scaleint/5 #because these slices are a fifth the size of the full total
    
        mapmap_scaled=mapmap/np.reshape(scaleint,mapmap.shape)
    
        mapmap_scaled=mapmap_scaled/lt_scaling[i-1]
    
    
        mapmap=mapmap_scaled
        mapmap_diff=mapmap_scaled
        
         
        np_lon2=north_plotting[0,:,3]
        np_lat2=north_plotting[0,:,2]
        np_loc2=north_plotting[0,:,1]
        
        
        
        
        fit_xmo,fit_ymo=fit_fieldline_fixed(np_lon2,np_lat2) 
        fity,fitx = fit_fieldline_scattered(ofl_lon,ofl_lat,sigma=6)
        
        
    
        
        mapmap2,lon2d2,lat2d2=masksubarray(lon2d,lat2d,mapmap,fit_xmo,fit_ymo)
        mapmap2_diff,lon2d2,lat2d2=masksubarray(lon2d,lat2d,mapmap_diff,fit_xmo,fit_ymo)
        
        
        
        mapmap3,lon2d3,lat2d3=masksubarray(lon2d2,lat2d2,mapmap2,fitx,fity,invert=True)
        mapmap3_diff,lon2d3,lat2d3=masksubarray(lon2d2,lat2d2,mapmap2_diff,fitx,fity,invert=True)
    
    
        #these dont actually change anything in the final plots
        # ellipse fit
        maglong=173.36
        maglat=71.54
        # by eye
        # maglong=175
        # maglat=72
        # Grodent 2003
        maglong=175
        maglat=75  
        
        
        
        lat2d4,lon2d4 =to_magnetic(lat2d3,lon2d3,maglong=maglong,maglat=maglat)
        fity2,fitx2=to_magnetic(fity,fitx,maglong=maglong,maglat=maglat)
    
    
    
    
        nplon1=north_plotting[3,:,3]
        nplat1=north_plotting[3,:,2]
        nplon2=north_plotting[13,:,3]
        nplat2=north_plotting[13,:,2]
        nplon3=north_plotting[23,:,3]
        nplat3=north_plotting[23,:,2]
        
    
        nplat1b,nplon1b=to_magnetic(nplat1,nplon1,maglong=maglong,maglat=maglat)
        nplat2b,nplon2b=to_magnetic(nplat2,nplon2,maglong=maglong,maglat=maglat)
        nplat3b,nplon3b=to_magnetic(nplat3,nplon3,maglong=maglong,maglat=maglat)
    
    
    
    # sector1_lat,sector1_lon,,transform=ccrs.PlateCarree(),c=sector1_loc
    
    
        sector1m_lat,sector1m_lon=to_magnetic(sector1_lat,sector1_lon,maglong=maglong,maglat=maglat)
        sector2m_lat,sector2m_lon=to_magnetic(sector2_lat,sector2_lon,maglong=maglong,maglat=maglat)
        sector3m_lat,sector3m_lon=to_magnetic(sector3_lat,sector3_lon,maglong=maglong,maglat=maglat)
        sector4m_lat,sector4m_lon=to_magnetic(sector4_lat,sector4_lon,maglong=maglong,maglat=maglat)
        sector5m_lat,sector5m_lon=to_magnetic(sector5_lat,sector5_lon,maglong=maglong,maglat=maglat)
    
    
    
        # if i==1:
        #     ax2.scatter(sector1m_lon,sector1m_lat,transform=ccrs.PlateCarree(),c=sector1_loc,vmin=0,vmax=24,cmap='hsv',marker='.')
        # if i==2:
        #     ax2.scatter(sector2m_lon,sector2m_lat,transform=ccrs.PlateCarree(),c=sector2_loc,vmin=0,vmax=24,cmap='hsv',marker='.')
        # if i==3:
        #     ax2.scatter(sector3m_lon,sector3m_lat,transform=ccrs.PlateCarree(),c=sector3_loc,vmin=0,vmax=24,cmap='hsv',marker='.')
        # if i==4:
        #     ax2.scatter(sector4m_lon,sector4m_lat,transform=ccrs.PlateCarree(),c=sector4_loc,vmin=0,vmax=24,cmap='hsv',marker='.')
        # if i==5:
        #     ax2.scatter(sector5m_lon,sector5m_lat,transform=ccrs.PlateCarree(),c=sector5_loc,vmin=0,vmax=24,cmap='hsv',marker='.')
    
    
    
        # plt.figure()
        # crs = ccrs.RotatedPole(globe=ccrs.Globe(flattening=(0.0)))
        # ax1 = plt.subplot(1,1,1,projection=ccrs.Orthographic(180, 80))
        # gl = ax1.gridlines(crs=ccrs.PlateCarree(), linewidth=1, color='grey', alpha=0.5, linestyle='dotted', draw_labels=True)
        # plt.tricontourf(lon2d2, lat2d2, mapmap3,levels=256,cmap='afmhot',vmax=maxint,vmin=minint) 
        # ax1.plot(fitx,fity,transform=ccrs.PlateCarree())
        # plt.show()
        
        
        # ax2.set_extent([0, 360, 50, 90], crs=ccrs.PlateCarree())
        
        # plt.scatter(sector2_lon,sector2_lat,transform=ccrs.PlateCarree(),c=sector2_loc,vmin=0,vmax=24,cmap='hsv',marker='.')
    
        # ax1.set_global()
        # 
        
        
        # ax2.plot(nplon1b,nplat1b,transform=ccrs.PlateCarree(),linewidth=1,alpha=0.5)
        # ax2.plot(nplon2b,nplat2b,transform=ccrs.PlateCarree(),linewidth=1,alpha=0.5,linestyle='--')
        # ax2.plot(nplon3b,nplat3b,transform=ccrs.PlateCarree(),linewidth=1,alpha=0.5,linestyle='--')
        
        
    # fig = plt.figure()
    # ax2 = plt.subplot(1,1,1,projection=ccrs.Orthographic(360-polelong, 90))
    
    
    
    
        if doallfive:
            if i==1:
        
                localtimes1=data_atlatlong(sector1_loc,sector1_lat,sector1_lon,lat2d3,lon2d3,smooth=True)    
                distance1=data_atlatlong(sector1_dis,sector1_lat,sector1_lon,lat2d3,lon2d3,smooth=True)    
                xsxs=np.isfinite(localtimes1)
                localtimes1=localtimes1[xsxs]
                distance1=distance1[xsxs]
                
                lat2d5=lat2d3[xsxs]
                lon2d5=lon2d3[xsxs]
                
                mapmap4_diff_loc=mapmap3_diff[xsxs]
                mapmap4_loc=mapmap3[xsxs]
                localtimes=localtimes1
                distance=distance1
        
        
            if i==2:
              
                localtimes1=data_atlatlong(sector2_loc,sector2_lat,sector2_lon,lat2d3,lon2d3,smooth=True)    
                distance1=data_atlatlong(sector2_dis,sector2_lat,sector2_lon,lat2d3,lon2d3,smooth=True)    
                xsxs=np.isfinite(localtimes1)
                localtimes1=localtimes1[xsxs]
                distance1=distance1[xsxs]
        
                lat2d5=np.append(lat2d5,lat2d3[xsxs])
                lon2d5=np.append(lon2d5,lon2d3[xsxs])
                localtimes=np.append(localtimes,localtimes1)
                distance=np.append(distance,distance1)
                mapmap4_diff_loc=np.append(mapmap4_diff_loc,mapmap3_diff[xsxs])
                mapmap4_loc=np.append(mapmap4_loc,mapmap3[xsxs])
        
            if i==3:
              
                localtimes1=data_atlatlong(sector3_loc,sector3_lat,sector3_lon,lat2d3,lon2d3,smooth=True)    
                distance1=data_atlatlong(sector3_dis,sector3_lat,sector3_lon,lat2d3,lon2d3,smooth=True)    
                xsxs=np.isfinite(localtimes1)
                localtimes1=localtimes1[xsxs]
                distance1=distance1[xsxs]
        
                lat2d5=np.append(lat2d5,lat2d3[xsxs])
                lon2d5=np.append(lon2d5,lon2d3[xsxs])
                localtimes=np.append(localtimes,localtimes1)
                distance=np.append(distance,distance1)
                mapmap4_diff_loc=np.append(mapmap4_diff_loc,mapmap3_diff[xsxs])
                mapmap4_loc=np.append(mapmap4_loc,mapmap3[xsxs])
        
            if i==4:
              
                localtimes1=data_atlatlong(sector4_loc,sector4_lat,sector4_lon,lat2d3,lon2d3,smooth=True)    
                distance1=data_atlatlong(sector4_dis,sector4_lat,sector4_lon,lat2d3,lon2d3,smooth=True)    
                xsxs=np.isfinite(localtimes1)
                localtimes1=localtimes1[xsxs]
                distance1=distance1[xsxs]
        
                lat2d5=np.append(lat2d5,lat2d3[xsxs])
                lon2d5=np.append(lon2d5,lon2d3[xsxs])
                localtimes=np.append(localtimes,localtimes1)
                distance=np.append(distance,distance1)
                mapmap4_diff_loc=np.append(mapmap4_diff_loc,mapmap3_diff[xsxs])
                mapmap4_loc=np.append(mapmap4_loc,mapmap3[xsxs])
        
            if i==5:
              
                localtimes1=data_atlatlong(sector5_loc,sector5_lat,sector5_lon,lat2d3,lon2d3,smooth=True)    
                distance1=data_atlatlong(sector5_dis,sector5_lat,sector5_lon,lat2d3,lon2d3,smooth=True)    
                xsxs=np.isfinite(localtimes1)
                localtimes1=localtimes1[xsxs]
                distance1=distance1[xsxs]
        
                lat2d5=np.append(lat2d5,lat2d3[xsxs])
                lon2d5=np.append(lon2d5,lon2d3[xsxs])
                localtimes=np.append(localtimes,localtimes1)
                distance=np.append(distance,distance1)
                mapmap4_diff_loc=np.append(mapmap4_diff_loc,mapmap3_diff[xsxs])
                mapmap4_loc=np.append(mapmap4_loc,mapmap3[xsxs])
    
     
        else:
            
            if i==1:
              
                localtimes1=data_atlatlong(sector1_loc,sector1_lat,sector1_lon,lat2d3,lon2d3,smooth=True)    
                distance1=data_atlatlong(sector1_dis,sector1_lat,sector1_lon,lat2d3,lon2d3,smooth=True)    
                
                xsxs=np.isfinite(localtimes1)
                localtimes1=localtimes1[xsxs]
                distance1=distance1[xsxs]
                
                lat2d5=lat2d3[xsxs]
                lon2d5=lon2d3[xsxs]
                
                mapmap4_diff_loc=mapmap3_diff[xsxs]
                mapmap4_loc=mapmap3[xsxs]
                localtimes=localtimes1
                distance=distance1
                
            
            if i==2:
              
                localtimes1=data_atlatlong(sector2_loc,sector2_lat,sector2_lon,lat2d3,lon2d3,smooth=True)    
                distance1=data_atlatlong(sector2_dis,sector2_lat,sector2_lon,lat2d3,lon2d3,smooth=True)    
                
                xsxs=np.isfinite(localtimes1)
                localtimes1=localtimes1[xsxs]
                distance1=distance1[xsxs]
                
                lat2d5=lat2d3[xsxs]
                lon2d5=lon2d3[xsxs]
                
                mapmap4_diff_loc=mapmap3_diff[xsxs]
                mapmap4_loc=mapmap3[xsxs]
                localtimes=localtimes1
                distance=distance1
                
            
            if i==3:
              
                localtimes1=data_atlatlong(sector3_loc,sector3_lat,sector3_lon,lat2d3,lon2d3,smooth=True)    
                distance1=data_atlatlong(sector3_dis,sector3_lat,sector3_lon,lat2d3,lon2d3,smooth=True)    
                
                xsxs=np.isfinite(localtimes1)
                localtimes1=localtimes1[xsxs]
                distance1=distance1[xsxs]
                
                lat2d5=lat2d3[xsxs]
                lon2d5=lon2d3[xsxs]
                
                mapmap4_diff_loc=mapmap3_diff[xsxs]
                mapmap4_loc=mapmap3[xsxs]
                localtimes=localtimes1
                distance=distance1
                
            
            if i==4:
              
                localtimes1=data_atlatlong(sector4_loc,sector4_lat,sector4_lon,lat2d3,lon2d3,smooth=True)    
                distance1=data_atlatlong(sector4_dis,sector4_lat,sector4_lon,lat2d3,lon2d3,smooth=True)    
                
                xsxs=np.isfinite(localtimes1)
                localtimes1=localtimes1[xsxs]
                distance1=distance1[xsxs]
                
                lat2d5=lat2d3[xsxs]
                lon2d5=lon2d3[xsxs]
                
                mapmap4_diff_loc=mapmap3_diff[xsxs]
                mapmap4_loc=mapmap3[xsxs]
                localtimes=localtimes1
                distance=distance1
                
            if i==5:
              
                localtimes1=data_atlatlong(sector5_loc,sector5_lat,sector5_lon,lat2d3,lon2d3,smooth=True)    
                distance1=data_atlatlong(sector5_dis,sector5_lat,sector5_lon,lat2d3,lon2d3,smooth=True)    
                
                xsxs=np.isfinite(localtimes1)
                localtimes1=localtimes1[xsxs]
                distance1=distance1[xsxs]
                
                lat2d5=lat2d3[xsxs]
                lon2d5=lon2d3[xsxs]
                
                mapmap4_diff_loc=mapmap3_diff[xsxs]
                mapmap4_loc=mapmap3[xsxs]
                localtimes=localtimes1
                distance=distance1
    
    
    
    
        if (not doallfive) | (i==5):
    
    
            if doallfive:
            # FULL SLICING HERE
                localtimes1=data_atlatlong(sector1_loc,sector1m_lat,sector1m_lon,lat2d4,lon2d4,smooth=True)
                localtimes2=data_atlatlong(sector2_loc,sector2m_lat,sector2m_lon,lat2d4,lon2d4,smooth=True)
                localtimes3=data_atlatlong(sector3_loc,sector3m_lat,sector3m_lon,lat2d4,lon2d4,smooth=True)
                localtimes4=data_atlatlong(sector4_loc,sector4m_lat,sector4m_lon,lat2d4,lon2d4,smooth=True)
                localtimes5=data_atlatlong(sector5_loc,sector5m_lat,sector5m_lon,lat2d4,lon2d4,smooth=True)
            
            
                for lts in range(36):
                    
                    slice_lat=np.array([90,60,30,0,-30,-60,-90,-60,-30,0,30,60,90])
                    slice_lon=np.array([lts*10,lts*10,lts*10,lts*10,lts*10,lts*10,lts*10,lts*10,(lts+1)*10,(lts+1)*10,(lts+1)*(lts+1)*10,(lts+1)*10,(lts+1)*10])
                    
                    
                    mapmap5_slice,lon2d5_slice,lat2d5_slice=masksubarray(lon2d4,lat2d4,mapmap3,slice_lon,slice_lat,invert=False)
                    mapmap5_diff_slice,lon2d5_slice,lat2d5_slice=masksubarray(lon2d4,lat2d4,mapmap3_diff,slice_lon,slice_lat,invert=False)
                    
                    
                    zaza=np.isfinite(mapmap5_diff_slice)
                    lon2d6_slice=lon2d5_slice[zaza]
                    lat2d6_slice=lat2d5_slice[zaza]
                    mapmap6_diff_slice=mapmap5_diff_slice[zaza]
                    
                    xvxv=np.nonzero(mapmap6_diff_slice)
                    lon2d7_slice=lon2d6_slice[xvxv]
                    lat2d7_slice=lat2d6_slice[xvxv]
                    mapmap7_diff_slice=mapmap6_diff_slice[xvxv]
             
                
            
                    # here, data is the local time values to be interpolated
                    # localtimes=data_atlatlong(data,datalat,datalon,newlat,newlon,smooth=True)
                    
            
                    xcvx=mapmap5_diff_slice[np.isfinite(mapmap5_diff_slice)]
                    xvxvx=xcvx[np.nonzero(xcvx)]
            
                    
                    mapmap5_diff2_slice=mapmap5_diff_slice
                    
                    xcvx2=mapmap5_diff2_slice[np.isfinite(mapmap5_diff2_slice)]
                    xvxvx2=xcvx2[np.nonzero(xcvx)]
                   
                    
                   # mapmap5_diff2mis the fully corrected intensities (except for a final scaling)
                   # just need to calculate local times for each lat/long of this data
                   # that is stored in an array called:sector1_loc, sector1_lat and sector1_lon
                   # what was that procedure for cross-correlating interp or something - look that up now....
                    # cant find it - was it on the PC?
                   
                    # print(i,lts,mapmap5_slice.size)    
                    
                    
                
                    # gl = ax2.gridlines(crs=ccrs.PlateCarree(), linewidth=1, color='grey', alpha=0.5, linestyle='dotted', draw_labels=True)
                
                    # ax2.scatter(lon2d5, lat2d5,transform=ccrs.PlateCarree(),c=mapmap5,cmap='afmhot',vmax=maxint,vmin=minint) 
                    # plt.show()
                    brightness[i,lts]= np.mean(mapmap5_slice)
                    brightness_diff[i,lts]= np.mean(xvxvx)
                    brightness_diff2[i,lts]= np.mean(xvxvx2)
            
        # fig.savefig('asd2.pdf', dpi=300, bbox_inches='tight', facecolor='white') 
        
        
        
            
                brightness_diff[np.isnan(brightness_diff)]=0
                
                
                # plt.show()
                # plt.imshow(brightness[1:5,:])
                b2=brightness+0.
                # print()
                b2_diff=brightness_diff+0.
                b2_diff2=brightness_diff2+0.
                
                b22=b2+0
                b2d2=b2_diff+0
                for y in range(5): b2[y+1,:]=b22[y+1,:]/np.mean(b22[1:5,:],axis=0)
                for y in range(5): b2_diff[y+1,5:30]=b2d2[y+1,5:30]/np.mean(b2d2[1:5,5:30],axis=0)
                b22=b2+0
                
                for y in range(5): b2[y+1,5:30]=b22[y+1,5:30]/np.mean(b22[y+1,5:30])
                
                b2[:,0:5]=0
                b2[:,30:]=0
                
                b2_diff[:,0:5]=0
                b2_diff[:,30:]=0
                
                b2_diff2[:,0:5]=0
                b2_diff2[:,30:]=0
                
                
                # plt.show()
                # plt.imshow(b2[1:5,:])
                # plt.show()
                # plt.imshow(b2_diff[1:5,:])
                # plt.show()
                # plt.figure()
                
                
                def zeromean(input):
                    step2=input[np.nonzero(input)]
                    # print(step2.size)
                    output = np.mean(step2)
                    return output
                
                xrange=360-np.arange(36)*10
                
                # for i in range(5):
                #     plt.plot(xrange,b2[i+1,:])
                #     # plt.ylim([0.6,1.4])
                
                # plt.figure()
                # for i in range(5):
                #     plt.plot(xrange,b2_diff[i+1,:])
                #     # plt.plot(b2_diff[i+1,:])
                #     # plt.ylim([0.6,1.4])
                
                
                
                # plt.figure()
                # for i in range(5):
                #     plt.plot(xrange,b2_diff2[i+1,:]/zeromean(b2_diff2[i+1,:]))
                #     # plt.ylim([0.6,1.4])
                
                
            
            
            
            
            
            
            
            
            
            nottoopolar=np.where(lat2d5 < 80)
            
            
            lt3=localtimes[nottoopolar]
            di3=distance[nottoopolar]
            mm3=mapmap4_diff_loc[nottoopolar]
            mm3=mapmap4_loc[nottoopolar]
            lon3=lon2d5[nottoopolar]
            lat3=lat2d5[nottoopolar]
            
            
            noheadwinds = np.where((lt3 < 6) | (lt3 > 18) | (di3 < 130) )
            
            # print(nottoopolar,noheadwinds)
            
            lt2=lt3[noheadwinds]
            di2=di3[noheadwinds]
            mm2=mm3[noheadwinds]
            # mm2b=mm3b[noheadwinds]
            lon2=lon3[noheadwinds]
            lat2=lat3[noheadwinds]
            
            
            # plt.figure()
            # plt.scatter(localtimes,distance,c=mapmap4_diff_loc,vmin=0.001,vmax=0.002,marker='.',lw=0.01)
            # plt.figure()
            # plt.scatter(lt2,di2,c=mm2,vmin=0.001,vmax=0.002,marker='.',lw=0.01)
            
            
            ltdeg=localtimes/24*360
            ltrad=ltdeg/360*np.pi*2
            # theta=ltrad* np.cos(ltrad)
            rad=distance
            
            
            
            # fig1, ax = plt.subplots()
            # ax.set_box_aspect(1)
            
            # plt.scatter(rad * np.cos(ltrad), rad * np.sin(ltrad),c=mapmap4_diff_loc,vmin=0.001,vmax=0.002,marker='.',lw=0.01)
            
            ltrad2=lt2/24*np.pi*2
            rad2=di2
            
            
            if doallfive:
                    pass
                
            else:
            
                if i == 1:
                    fig, axs = plt.subplots(2, 5, gridspec_kw={'width_ratios': [1,1,1,1,1]}) 
              
                    fig.set_figwidth(15*1.5)
                    fig.set_figheight(3*3)
                    fig.set_dpi(100)
                    fig.subplots_adjust(wspace=0, hspace=0)
                
                # axs[0,5-i].set_box_aspect(300/250)
                # axs[0,5-i].scatter(rad2 * np.cos(ltrad2), rad2 * np.sin(ltrad2),c=mm2,marker='.',lw=0.01,vmin=0.4,vmax=1.4)
                
                
                # print('max mm = ',np.max(mm2))
                
                # plot_labels_on(axs[0,5-i],white=False,inum=i)
        
                
                axs[0,5-i].set_box_aspect(300/250)

                
                cc=axs[0,5-i].tricontourf(rad2 * np.cos(ltrad2), rad2 * np.sin(ltrad2),mm2, levels=np.linspace(0.4,1.2,128),cmap='inferno')
                axs[0,5-i].tricontour(rad2 * np.cos(ltrad2), rad2 * np.sin(ltrad2),mm2, levels=np.linspace(0.4,1.2,128),linewidths=1,cmap='inferno')
                # cc=axs[0,5-i].scatter(rad2 * np.cos(ltrad2), rad2 * np.sin(ltrad2),c=mm2)
                axs[0,5-i].set_xlim([-100,150])
                axs[0,5-i].set_ylim([-150,150])
                
                axs[0,5-i].text(-70,-130, lt_text[5-i] , ha="center", va="center",c='k',zorder=4,fontsize='large',bbox=bbox_props1)

                # if (not doallfive):
                #     ilabel="PLT: "+str(18-i*2)
                # else:
                ilabel=""
                plot_labels_on(axs[0,5-i],ilabel=ilabel,inum=i)
                
                
                if i == 5:
                    axs[1,0].axis('off')
                    # axs[1,1].axis('off')
                    # axs[1,2].axis('off')
                    # axs[1,3].axis('off')
                    axs[1,4].axis('off')
                    ax_cb=plt.subplot(5,3,8+3)
                    ax_cb.axis('off')

                    cbar=fig.colorbar(cc,ax=ax_cb,location='top', ticks=[0.4,0.6,0.8,1.0,1.2],aspect=40)
                    cbar.ax.set_xlabel('Auroral intensity')
                    cbar.ax.xaxis.set_label_position("bottom")
                    # # axs[5].set_yticks([0.4,0.6,0.8,1.0,1.2,1.4])
                    # axs[5].set_ylim([0.4,1.4])
                    fig.savefig('asd_lt_di_indi.pdf', dpi=300, bbox_inches='tight', facecolor='white') 

                    
                # txt4.set_bbox(dict(facecolor='k', alpha=0.3, edgecolor='k', boxstyle='round,pad=0.2'))
                
                
                # fig1, ax = plt.subplots()
                # ax.set_box_aspect(1)
                
            
            # plt.tricontourf(rad2 * np.cos(ltrad2), rad2 * np.sin(ltrad2),mm2b,vmin=0.001,vmax=1.7,levels=256)
            
            
            XJ=rad2 * np.cos(ltrad2)
            YJ=rad2 * np.sin(ltrad2)
            
            xx=rad2 * np.cos(ltrad2)
            yy=rad2 * np.sin(ltrad2)
            
            # from matplotlib.patches import Circle
            # x0=0
            # y0=0
            # radii=15
            # patches = []
            # circle = Circle((x0, y0), radii,color='w')
            # ax.add_patch(circle)
            
            
            
            
            mm5=mm2+0
            mm5[:]=0
            
            mm4=(mm2*-1)
            mm4=mm4**10
            # mm4[:]=0
            
            mm5[np.sqrt(XJ**2+YJ**2) < 30] = 240
            for i in range(0,40,2): mm5[(lon2 >140+i) & (lon2<142+i) & (lat2 > 70-i/4) & (lat2 < 75-i/4)] = 120
            
            
            # print(XJ[mm4 > np.max(mm4/2)],YJ[mm4 > np.max(mm4/2)])
            
            
            
            # mm4[mm2 < np.max(mm2)/10] = 240
            
            
            
            
            xx=rad2 * np.cos(ltrad2)
            yy=rad2 * np.sin(ltrad2)
            
            gridded_2dint=np.zeros([60,50])
            gridded_XX=np.zeros([60,50])
            gridded_YY=np.zeros([60,50])
            gridded_2dstdev=np.zeros([60,50])
            gridded_2dcount=np.zeros([60,50])
            
            for xxx in range(50):
                for yyy in range(60):
            
                    xmin=xxx*5-100
                    xmax=(xxx+1)*5-100
                    ymin=yyy*5-150
                    ymax=(yyy+1)*5-150
            
                    truevalues=mm2[(xx >= xmin) & (xx <= xmax) & (yy >= ymin)& (yy <= ymax)]
                    truexx=xx[(xx >= xmin) & (xx <= xmax) & (yy >= ymin)& (yy <= ymax)]
                    trueyy=yy[(xx >= xmin) & (xx <= xmax) & (yy >= ymin)& (yy <= ymax)]
            
                    # print(xxx,yyy,np.size(truevalues))
            
                    if np.size(truevalues) > 0:
                        gridded_2dint[yyy,xxx]=np.median(truevalues)
                        gridded_XX[yyy,xxx]=np.mean(truexx)
                        gridded_YY[yyy,xxx]=np.mean(trueyy)
                        gridded_2dstdev[yyy,xxx]=np.std(truevalues)
                        gridded_2dcount[yyy,xxx]=truevalues.size
                        
                      
            x = np.arange(-100, 150, 50)
            y = np.arange(-150, 150, 60)
            X, Y = np.meshgrid(x, y)
     
            
    
    
    if doallfive:
        # fig2, axs2 = plt.subplots(1,5)
        fig2, axs2 = plt.subplots(1,3)
    
    
        fig2.subplots_adjust(wspace=.2, hspace=5)

        fig2.set_figwidth(12)
        fig2.set_figheight(6)
        fig2.set_dpi(100)
     
        # axs2[0].set_box_aspect(300/250)
    
        # axs2[0].scatter(rad2 * np.cos(ltrad2), rad2 * np.sin(ltrad2),c=mm2,marker='.',lw=0.01,vmin=0.4,vmax=1.4)
        # axs2[0].hexbin(rad2 * np.cos(ltrad2), rad2 * np.sin(ltrad2),vmin=0.4,vmax=1.4,gridsize=10,mincnt=5)
        # axs2[0].hexbin(rad2 * np.cos(ltrad2), rad2 * np.sin(ltrad2),vmin=0.4,vmax=1.4,gridsize=20,mincnt=2)
        # axs2[0].hexbin(rad2 * np.cos(ltrad2), rad2 * np.sin(ltrad2),vmin=0.4,vmax=1.4,gridsize=40,mincnt=2)
        # axs2[0].hexbin(rad2 * np.cos(ltrad2), rad2 * np.sin(ltrad2),vmin=0.4,vmax=1.4,gridsize=80,mincnt=2)
        
        
        print('max mm = ',np.max(mm2))
        
    
        
        
        axs2[0].set_box_aspect(300/250)
    
        # axs2[0].tricontourf(rad2 * np.cos(ltrad2), rad2 * np.sin(ltrad2),mm2, levels=np.linspace(0.4,1.4,64))
        # axs2[0].tricontour(rad2 * np.cos(ltrad2), rad2 * np.sin(ltrad2),mm2, levels=np.linspace(0.4,1.4,64),linewidths=1)
        # plot_labels_on(axs2[0],white=False)

    
        axs2[1].set_box_aspect(300/250)
        axs2[2].set_box_aspect(300/250)
        # axs2[4].set_box_aspect(300/250)
    
    
    
        trigrid_2dint=gridded_2dint[gridded_2dint>0]
        trigrid_2dstdev=gridded_2dstdev[gridded_2dint>0]
        trigrid_2dcount=gridded_2dcount[gridded_2dint>0]
        trigrid_XX=gridded_XX[gridded_2dint>0]
        trigrid_YY=gridded_YY[gridded_2dint>0]
    
        # axs2[2].imshow(gridded_2dint,origin='lower', extent=(-100, 150, -150, 150)) 
        # axs2[3].imshow(gridded_2dstdev,origin='lower', extent=(-100, 150, -150, 150))  
        # axs2[4].imshow(gridded_2dcount,origin='lower', extent=(-100, 150, -150, 150),vmax=10)  
          
    
        cc1=axs2[0].tricontourf(trigrid_XX, trigrid_YY,trigrid_2dint, levels=np.linspace(0.4,1.2,128),cmap='inferno')
        axs2[0].tricontour(trigrid_XX, trigrid_YY,trigrid_2dint, levels=np.linspace(0.4,1.2,128),linewidths=1,cmap='inferno')
        plot_labels_on(axs2[0],white=True)

    
        cc2=axs2[1].tricontourf(trigrid_XX, trigrid_YY,trigrid_2dstdev, levels=np.linspace(0.0,1.5*0.2,128),cmap='jet')
        axs2[1].tricontourf(trigrid_XX, trigrid_YY,trigrid_2dstdev, levels=np.linspace(0.0,1.5*0.2,128),linewidths=1,cmap='jet')
        plot_labels_on(axs2[1],white=True)
        cc3=axs2[2].tricontourf(trigrid_XX, trigrid_YY,trigrid_2dcount, levels=np.linspace(0,20,128),cmap='gnuplot')
        axs2[2].tricontourf(trigrid_XX, trigrid_YY,trigrid_2dcount, levels=np.linspace(0,20,128),linewidths=1,cmap='gnuplot')
        plot_labels_on(axs2[2],white=True)
    
        # axs2[0].axis('off')
        cbar=fig2.colorbar(cc1,ax=axs2[0], ticks=[0.4,0.6,0.8,1.0,1.2],location="bottom", pad=0.1)
        cbar.ax.set_xlabel('Auroral intensity')
        cbar.ax.xaxis.set_label_position("top")

        # axs2[1].axis('off')
        cbar=fig2.colorbar(cc2,ax=axs2[1], ticks=[0,0.1,0.2,0.3],location="bottom", pad=0.1)
        cbar.ax.set_xlabel('Std. dev of emission')
        cbar.ax.xaxis.set_label_position("top")

        # axs2[2].axis('off')
        cbar=fig2.colorbar(cc3,ax=axs2[2], ticks=[0,5,10,15,20],location="bottom", pad=0.1)
        cbar.ax.set_xlabel('Counts per 10R$_J$$^2$')
        cbar.ax.xaxis.set_label_position("top")
   
    
    
        fig2.savefig('asd_lt_di_total.pdf', dpi=300, bbox_inches='tight', facecolor='white', pad_inches=0) 
      
                    
